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Abstract 

The bifurcation diagram of a model stochastic differential equation with delayed feedback is pre- 
sented. We are motivated by recent research on stochastic effects in models of transcriptional gene 
regulation. We start from the normal form for a pitchfork bifurcation, and add multiplicative or 
parametric noise and linear delayed feedback. The latter is sufficient to originate a Hopf bifurcation 
in that region of parameters in which there is a sufficiently strong negative feedback. We find a 
sharp bifurcation in parameter space, and define the threshold as the point in which the stationary 
distribution function p{x) changes from a delta function at the trivial state x = to p{x) ~ 
at small x (with a = — 1 exactly at threshold). We find that the bifurcation threshold is shifted 
by fluctuations relative to the deterministic limit by an amount that scales linearly with the noise 
intensity. Analytic calculations of the bifurcation threshold are also presented in the limit of small 
delay r ^ that compare quite favorably with the numerical solutions even for r = 1. 

PACS numbers: 02.30.Ks, 05.10.Gg, 05.70.Ln, 87.16.Yc, 87.18. Cf 



I. INTRODUCTION 



We study the bifurcation diagram of a nonlinear stochastic differential equation that 
includes delayed feedback. The model equation considered exhibits both pitchfork and Hopf 
bifurcations. The bifurcation thresholds are obtained as a function of the model parameters, 
and our results contrasted with two related limits: The deterministic limit of a differential 
delay equation, and the stochastic bifurcation of the same model without delay. 

The study of differential delay equations [ll] is an important topic in applied mathe- 
matics, with widespread applications in Physics (lasers, liquid crystals), control systems in 



Physiology (neural and cardiac tissue activity) 



and Economy (agricultural commodity 



prices) j^. Recent interest has arisen in the mathematical modeling of cellular function at 
the molecular level, especially in transcriptional gene regulation ^. Feedback regulation is 
a common motif in cellular networks, with delays arising from the complexity of the un- 
derlying network, or from the wide disparity in time scales of the many chemical processes 
involved in regulation For example, DNA is transcribed at a rate of 10 to 100 nucleotides 
per second, and it may take a delay of the order of minutes before the transcription factor 
appears as a finished product in the cell, and hence is available for regulation. Significant 
delays can also be attributable to the time required for the diffusion of proteins through 
membranes, so that, for example, the auto regulated feedback on protein production at time 
t is often proportional to protein concentration at time t — r, where r is known as the delay 
time. For short delay times, a reaction may be approximated as being instantaneous, and 
the system as being in quasi equilibrium. However, when the delay is comparable to the 
characteristic time scales of reaction, the non instantaneous nature of the interactions can 



no longer be ignored, and de 



the network under study 



a. 



ay terms need to be included in the governing equations for 



Experimental evidence has been mounting that highlights the importance of stochastic 



effects in transcriptional regulation [?, 8 
neered gene circuits and networks as well 



9. 



IG 
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, not only for natural networks, but for engi- 



12l |. However, despite the wealth of evidence 



pointing to the importance of stochasticity in feedback regulation, delays in stochastic mod- 
els of metabolic feedback are very often neglected, possibly because the resulting stochastic 
equations are no longer Markovian, and hence rarely tractable analytically. Exceptions in- 
clude the derivation of a two time Fokker-Planck equation and the study of its small delay 



2 



time limit in 13| , and results on the bifurcation of the first and second moments of a linear 
stochastic equation with delay [g], 1^. We extend these latter results to the analysis of the 
stationary probability distribution function of a nonlinear model, and discuss in detail the 
stability of the solutions that results from the interplay of delay and stochasticity. 

We focus here on the normal form for a pitchfork bifurcation augmented with multiplica- 
tive or parametric noise (additive noise has been shown to have no effect in the case of a 

n 

delayed equation [6|), and linear delayed feedback. The latter is sufficient to originate a 
Hopf bifurcation in some region of parameters (sufficiently strong negative feedback). Our 
model for the dynamical variable x{t) is 

x{t) = ax{t) + hx{t - r) - x^{t) + i{t)x{t) , (1) 

where the constant a plays the role of the control parameter, h is the intensity of a feedback 
loop of delay r, and ^(t) is a white, Gaussian noise of intensity D. The initial condition is 
a function (pit) specified on t = [— r, 0]. We generally focus on the stationary probability 
distribution function p{x) for a range of values of a, h and D. 

The bifurcation diagram for the differential delay equation resulting from Eq. ([T]) with 
^ = is known (see, e.g. [6,]). Linearization around x = shows that trajectories decay 
asymptotically to zero according to 

h< -a if 6r > -1 , (2) 



rV62 -a2 < cos-H-f) if &r < -1 . (3) 

The boundary separating exponentially decaying solutions from exponentially growing so- 
lutions is shown as the solid line in Fig. (jlj). The upper branch corresponds to a pitchfork 
bifurcation (real eigenvalue, Eq. ([21)), whereas the lower branch corresponds to a Hopf bi- 
furcation (complex eigenvalue, Eq. ([3])). In both cases, we show in the figure the case of 
r = 1. The cusp at the intersection of both boundaries is located at (a, 6) = (1/r, — 1/r). 

The bifurcation threshold of Eq. (fll) without delay ih = 0) is also known, and has been 
given in [15|, [161, [ITJ • Recall that an analysis of the linearized equation leads to the unphysical 
conclusion that the bifurcation threshold depends on the order of the statistical moment of 
x{t) considered. On the other hand, with a saturating nonlinearity in Eq. ([T]), a stationary 
probability distribution of x exists both below and above threshold, thus allowing a proper 
determination of the bifurcation point. The stationary distribution function of x with 6 = 



has been found to be [1 



a < —1 po{x) = 6{x) , (4) 
a > —1 Po{x) = Nx^e ^ (5) 

where the exponent a = a/ D — 1, and is a normahzation constant. The solution ([5]) exists 
but is not normahzable for a < —1 and hence it is not a physically admissible solution. 
Therefore the stochastic bifurcation threshold is located at = 0, point at which p{x) 
changes from a delta function at the the origin to a power law at small x with an exponential 
cut off at large x. Contrary to what an analysis of the moments of x{t) would indicate, the 
existence of parametric fluctuations has no effect on the location of the bifurcation threshold: 
both deterministic and stochastic equations exhibit a pitchfork bifurcation at Oc = 0. We 
further note that in — 1 < « < 0, p{x) is unimodal with a singularity at x = 0, whereas for 
a > the distribution is bimodal reflecting nonlinear saturation of x. 



II. STOCHASTIC BIFURCATION 



We now turn to the case of delay, b ^ 0. Analytical results for the stability of the trivial 
solution a; = for the linearization of Eq. ([1]) have been given in [g], and are shown 
in Fig. (jlj). The bifurcation threshold of the first moment (x) is shifted relative to the 
deterministic delay equation result; only bounds can be given for the stability of the second 
moment (x^) 1^, and no results are available for p{x). Given the anomalous behavior 
described above for the stochastic bifurcation of the linear equation with 6 = 0, it is of 
interest to determine p{x) for the full model of Eq. ([1]). Unfortunately, the non Markovian 
character of this equation has precluded progress along these lines 13 1. 

We have first extended an existing second order algorithm for the integration of stochastic 
differential equations 3] to the case of delay. The algorithm needs to take into account 
trajectories into the past for an interval r, and also new contributions from the stochastic 
terms that result from the coupling to the delayed feedback. Derivation of the algorithm is 
presented in appendix [Al In the numerical calculations to be presented below, the initial 
condition is a constant function in [— r, 0] for each trajectory, with the constant being drawn 
from a Gaussian distribution of zero mean and variance 1. The time step used in the 
numerical integration is At = 0.01. We also present approximate analytic calculations in 
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the limit of small delay time r ^ 0, following the approach of Frank 19(]. Reasonable 
agreement is found with our numerical calculations for r = 1. 



A. Pitchfork bifurcation 

A qualitative view of the bifurcation of Eq. ([1]) is given in Fig. ([1]). Equation ([1]) has 
been integrated numerically, and histograms of x{t) computed once trajectories approach a 
statistical steady state. The histograms shown correspond to 10^ independent trajectories 
for each value of a. For a ^ — 1, the histogram is approximately a delta function at a; = 0. 
As discussed further below, we observe long transients in x{t) until it eventually decays to 
X = 0. At a critical value Oc ~ —1, the bifurcation point, a broad distribution emerges, 
although the most likely value remains a; = 0. At larger values of a, the histogram becomes 
bimodal. The histogram shown in the figure corresponds to the direct bifurcation branch, 
but a qualitatively similar graph is obtained for the Hopf bifurcation. 

Our results for the stationary distribution function p{x) are shown in Fig. ([2]). For 
a < ac, we expect p{x) = 6{x). We observe instead a very long lasting transient with 
p{x) approximately characterized by a power law distribution, with an effective exponent 
a < — 1 (leading to a non normalizable distribution). The amplitude of the point p{x ^ 0) 
(not shown in the figure) grows with time, signaling the build up of the delta function at 
the origin. Because of overall normalization, growth at x = implies a decaying amplitude 
for finite x, as shown in the figure. For a > ac, we do obtain a time independent power 
law distribution p{x) with exponent — 1 < a < 0. This distribution is normalizable, and 
represents the stationary distribution above threshold. We finally show p{x) in the range of 
a for which the distribution is bimodal. The function p{x) around the most likely value is 
approximately constant over time, but we still observe some transients in the region around 
X = 0. Figure ([3]) shows our results for the exponent of a power law fit to p{x) at small x, 
as a function of the value of the control parameter a. We observe a smooth variation of the 
exponent a with a that allows a convenient determination of a^, the value of a for which 
a = —1. This is the method that we have used to determine the bifurcation threshold in all 
the results presented below. 

We summarize our results for the bifurcation diagram of Eq. ([T]) in Fig. (j3]). The 
analytic results for the threshold without noise (^ = 0) are shown for reference, as well as 
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earlier results for the threshold of (x) of the linearized equation 61. Our numerical results 

n 

do agree with the known threshold for the special point of no delay 6 = given in |15| . 
The figure presents our results for the stochastic bifurcation threshold defined directly from 
the stationary probability distribution function as discussed above. We conclude that the 
stochastic threshold is shifted relative to the deterministic threshold except in the special 
point of no delay (6 = 0). We have also examined in Fig. ([5]) the dependence of the threshold 
shift as a function of the noise intensity D. In analogy with the case of no delay, we find a 
linear dependence of Oc with D. 



B. Hopf bifurcation 

The calculation of p{x) just shown for the case of a pitchfork bifurcation has been repeated 
in the vicinity of the deterministic Hopf branch shown in Fig. (jlj). In the deterministic case, 
the bifurcation leads to oscillation. When fluctuations are added, oscillation amplitudes 
fluctuate as well. We also observe in this range of parameters a sharp bifurcation threshold 
which is shifted relative to the deterministic Hopf bifurcation. The bifurcation threshold is 
verified with the probability distribution function of the maximum amplitude of the Fourier 
transform of the trajectories. For each trajectory, we calculate the Fourier transform of 
the stationary solution over a finite window in time, identify its maximum amplitude, and 
construct an histogram over those maxima along the trajectory and over the ensemble. 

Our results for the stationary probability distribution function of the maximum amplitude 
of the Fourier transform are shown in Fig. ([6]) for three increasing values of the control 
parameter a. For a < Oc we observe an effective power law with exponent a < —1 and 
hence not a physical distribution. As was the case above, this is manifested by a transient 
distribution that decreases with time. As the value of a is increased, a power law distribution 
is found with exponent — 1 < a < 0. The distribution obtained is integrable and stationary. 
For yet larger values of a, the distribution becomes bimodal. 

Figure ([7]) shows the results of a power law fit to the resulting distributions at small x 
for a range of values of a of the two probability distribution functions introduced. We have 
undertaken this analysis for a range of values of b, and the resulting Hopf branch of the 
bifurcation diagram is shown in Fig. (jl]). Note that it is also shifted relative to the branch 
corresponding to the deterministic equation. 
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III. FOKKER-PLANCK EQUATION 



We next turn to an approximate analytic calculation of the stationary distribution func- 
tion p{x) for Eq. ([T]). The difficulty in obtaining a closed, analytic expression lies in 
the non Markovian nature of Eq. ([T]) and the associated need to find the joint proba- 
bility distribution p{x{t),x{t — r)). When r is larger than the correlation time of x, the 
derivation is simplified by assuming statistical independence between x{t) and x{t — r), or 
p{x{t),x{t — r)) = p{x{t))p{x{t — r)). This approximation has already been considered in 
the literature (e.g., ref. ^]). However, the assumption of independence does not hold near 
a bifurcation since the characteristic correlation time diverges |20|. We instead proceed as 
follows: Define the probability distribution of x as p{x,t) = {6{x{t) —x)). By using the 
properties of the Dirac delta function, and Eq. ([1]) one finds 
d d 

—p(x,t) = -—\(S(x(t)-x)(ax + bXr-x^)) + {S(x(t)-x)x^(t))] , (6) 
ot ox 

where we have introduced the dummy variables x{t) x and x{t — r) ^ Xr- The second 
term in the right hand side of Eq. ([6]) can be written as 

{6{x{t) - x)xm = ^(t) {6{x{t) - x)e(t)) , (7) 

with 

by using the Furutsu-Novikov theorem. By using the properties of the delta function, one 
can write 

and since = x{t) from Eq. ([1]), we find 

{5{x{t) - x)at)) = -D^xit)p{x, t) , (10) 

where we have also used 

/oo 
5{x{t) - x)f{x)p{x, t)dx = f{x{t))p{x, t) . (11) 
-oo 

Thus the second term in the right had side of Eq. is 



-^mxit)-x)xmm = D-^ 



d 

xit)—x{t)p{x,t) 



(12) 



[x{t)p{x, t)] + [x\t)p{x, t)] . 



Consider now the first term in the right hand side of Eq. ([6]). By using the law of total 
expectation, 



((5(x(t) — x){ax + bxr — x^)) 




5{x{t) — x){ax + bXr — x^)p{x, Xr)dxdxr , (13) 



We anticipate that independence, p{x,Xr) = p{x)p{xr), does not hold near the bifurca- 
tion threshold. We write instead p{x,Xr) = p{xr\x)p{x), where p{xr\x) is the conditional 
probability of finding Xr aX t — r given the information that x{t) = x. Then 



(^6{x{t) — x){ax + bXr — x^)) ~ J J ^{^{t) ~ ^){(^^ + bXr — x^)p{xr\x)p{x)dxdxr 

= p{x,t) / [ax{t) + bx-r — x^(t)] p{xr\x{t))dxT- 



ax{t) — x^{t) + b Xrp{xr\x{t))dxT- 



p{x,t) 



The last integral represents the conditional expected value of x{t — r) given x(t), 
Substitution of Eqs. (fT2l) and (IMI) into Eq. ([6]) leads to the Fokker-Planck equation. 



(14) 



(15) 



^p(a;, t) = [{a + D)x{t) - x\t) + b{xr\x{t))] p{x, t) + D-^ [x^it)pix, t)] . (16) 

We have not been able to determine {xr\x{t)) analytically. However, it is possible to 
derive an approximate expression under the assumption that the delay r is small 19|. We 
write drift terms in Eq. ([T]) in the Ito representation (recall to we were working with the 
Stratonovich representation), f{x,Xr) = {a+D)x—x^+bXr [2l|. Define /(x) = {a+D)x—x^, 
/(o)(x) = f{x) + Bx, and g{x) = X, so that f{x, Xr) = f{x) + Bxr- Furthermore, since we are 
in the small time delay approximation, one can assume that the zero-th order approximation 
pfj{xr,t — T\x,t) of the stationary conditional probability distribution function pst{xT,t — 
r|x,t) is Gaussian j22| and given by, 

(0), I I ^ K-x-/W(x)r]^^ 

Pst l-^T) ^ '-j 



1 



exp 



2TTTg'^{x) \ 2Tg^{x 
The approximate drift term of the Fokker-Planck equation is then, 

[Xr 



(17) 



fix) 



1 

2'KTX'^ 



dXrf{x, Xr) exp 



X 



/W(x)r]- 



(18) 
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Integrating, one finds, 

f{x) = (1 + 6r)/(°) (a;) = (1 + br) [{a + b + D)x - x^] . (19) 

An approximate expression for the conditional average of Xr given x{t) is obtained by com- 
paring Eq. ([in]) with Eq. ([H]), 

{xr\x{t)) = [1 + r(a + 6 + D)] x{t) - rx^ . (20) 

Further substitution of Eq. (120|) into Eq. (fT6l) results in a closed form for the steady state 
distribution, 

^) = - (TTfe^^ [{{a + b + D) x{t) - x'{t) } p{x, t)] + ^-^-^ A_ [x\t)p{x, t)] , 

(21) 

with stationary solution, 

p(a:) = A^|a;| d 20 ^ . (22) 

This stationary distribution is found to agree quite well with our numerical determination of 
p{x) (Fig. ([2])) around the pitchfork branch, but fails around the Hopf branch. Furthermore, 
we can now introduce an analytic determination of the bifurcation threshold as the point 
in which p(x) (Eq. ( l22l) ) ceases be normalizable: Oc = ""+^[^+^1^'^+^+-^)] _ 1 = _i. xhe 
bifurcation threshold is located at 

This prediction is also in very good agreement with our numerical results for the pitchfork 
bifurcation branch. 

We have also verified our results for the conditional average with a direct numerical 
integration of the model equation. We show our results for {xr\x{t)) as a function of x{t) 
in Fig. ([8]). Equation ([1]) is integrated numerically until it reaches a statistical stationary 
state. For every value of the dynamical variable x, the average of its value at time t — r is 
collected for 1 x 10^ independent trajectories in a time window t G (290, 300). We note that 
for large values of x, the results become less accurate because of fewer data points in this 
region. This analysis has been repeated for a range values of a, b, and D. As shown in the 
figure, we observe good agreement between the numerical results and Eq. (!20l) in the region 
around a; = 0. 
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In summary, we find that the bifurcation point in our stochastic differential equation 
with delay remains sharp. This is not a straightforward observation since the delay term 
in Eq. ([T]) could effectively act as an additive source of noise, and lead to an imperfect 
bifurcation instead. This does not appear to be the case. We also note that long transients 
can be expected below the bifurcation threshold, as all trajectories eventually decay to the 
trivial solution x = 0. Second, we observe that different moments of x obtained from the 
linearized equation bifurcate at different values of the control parameter. When a saturating 
nonlinearity is introduced into the model, all moments bifurcate at Oc. By defining the 
bifurcation point in the stochastic case as that in which the power law form of p{x) becomes 
normalizable (a = — 1), we show that the bifurcation threshold is shifted relative to that 
of the underlying deterministic equations. The shift goes to zero as 6 — > 0, and otherwise it 
scales linearly with the noise intensity D. We have also derived and approximate expression 
for the stationary distribution function p{x) in the limit of small delay time r. The threshold 
location obtained from this approximation agrees well with our numerical determination even 
when r = 1. 
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APPENDIX A: ALGORITHM FOR STOCHASTIC DIFFERENTIAL EQUA- 
TIONS WITH DELAY 

We summarize in this appendix the numerical algorithm developed to integrate Eq. ([T]). 
Formal integration yields, 

pt+At pt' pt+At pt+At 

x{t + At) = x{t) + ax{t')dt' - x^{t')dt'+ x{t')^{t')dt' + hx{t'-T)dt'. 

Jt Jt Jt Jt 

(Al) 

with 

x{t') = x{t) + J' ax{t")dt" - x\t")dt" + j'^ x{t")i{t")dt" + j' bx{t" - T)dt". (A2) 

We now approximate x{t") ~ x{t) and x{t" — r) ^ x{t — r) so that to first order in {t' — t) 
we have, 

x{t') = x{t) + [ax{t) + bx{t - r) - x^t)] {f -t) + x{t) i{t")dt" . (A3) 

The nonlinear term x'^{t') also must be expanded around x{t) 

x\t') ^ x{tf + 3a;2(t)(x(t') - x{t)) . (A4) 

We substitute Eqs. (]A3|) and (1A4|) into equation Eq. (lAll) and find, 

J.2 rt+At ft' 



{t + At) - x{t) = ax{t)At + [a'^x{t) + abx{t - r) - ax^{t)] h ax{t) / / ^{t")dt"dt' 

2 Jt Jt 

/t+At t-t+At 
i{t')dt' + [ax{t) + hx{t - r) - x^(t)] {t' - t)^it')dt' 

/t+At pt' 
i{f) ^{t")dt"dt' + bx{t-T)At 



+ x{t) 



Af 

+ [bax{t - r) + b^x{t - 2r) - bx^{t - r)] — 

pt+At pt'-T 

+ bx{t-T) j j i{t")dt"dt' -X^{t)At 

Jt Jt-T 

A .2 pt+At pt' 

- 3x\t) [ax{t) + bx{t - r) - x^{t)] 3x^{t) / / ^{t")dt"dt'. 

2 Jt Jt 

In order to calculate the integrals containing the random process ^(t), we define 

ft+At 

i{t')dt' = Gi{t,At) , (A5) 
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/ / i{t")dt"dt' = G2{t,At) . (A6) 
Jt Jt 

If ^{t) is a Gaussian process of zero mean, Gi and G2 are also Gaussian variables of zero 
mean, and correlations 

(Gl) = 2DAt , (AT) 
(Gl) = ^At3 , (A8) 
{G1G2) = DAt^ . (A9) 
The three remaining integrals can be expressed in terms of Gi and G2 as 

{t' - t)^{t')dt' = Gi{t, At) At - G2{t, At) , (AlO) 

I m I i{t")dt"dt' = Ai))^ , (All) 

/ / i{t")dt"dt' ^G2{t- T, At) . (A12) 

Jt Jt-T 

The Gaussian variables Gi and G2 can be simulated with two Gaussian random variables, 
^'i and ^2, of zero mean and variance one: 



Gi{t, At) = V2DAt^i{t) , (A13) 



G2{t, At) = ^^At3 i^^Mt) + lMt)j ■ (A14) 
Combining all our results, we write the iteration of our algorithm, 

x{t + At) = x{t) [l + aAt + c?— + (1 + aAt)Gi{t, At) + -{Gi{t, At))' 

( At^ \ 

+ x{t - t) [ bAt + '^ab— + bAtGi{t, At) - bG2{t, At) + bG2{t - r. At) \ 

+ x{t - 2t) i b^— j + x^{t) {-2aAt^ - {Gi{t, At) + l)At - 2G2) 

2,,, ,f?,bAe\ 3, ,fbAt^\ 5,,/3At2\ 
- x\t)x(t - t) y-^) - - -r) [-^j + x\t) j . (A15) 
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FIG. 1: Long time histogram of x (in grey scale) as a function of the control parameter a with 
6 = 1, r = 1, and D = 0.3. The histograms have been collected in the time interval t € (50,80) 
and further averaged over 1 x 10^ independent runs. In the absence of noise, the critical value of 
the control parameter for instability is Oc = — 1. We find instead that the bifurcation from a delta 
function to a power law distribution occurs at Oc =^ —1.18 instead for this set of parameters. 
Fluctuations around a; = are observed for a <ac due to the finite length of the time series. 
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FIG. 2: Probability distribution function p{x) for 6 = 1, r = 1 and D = 0.3 at the times given, 
and averaged over 10^ independent realizations. Values of the control parameter shown are: (a) 
a = -1.2 with a ~ -1.21, (b) a = -1.1 with a ~ -0.66, and (c) a = -0.9 with a ~ 0.61. The 
distributions in (a) show a clear transient, whereas those in (b) and (c) are stationary. The solid 
line shows the power law at small x; the domain covered by the line indicated the range of data 
that were used to estimate a, and is placed above or below the curves for clarity. The dashed line 
is our approximate determination of in the limit of small r. 
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FIG. 3: Results of a power law fit to p{x) for small x. The fitted value of the exponent a is shown 
as a function of the control parameter a. We define the bifurcation threshold for the stochastic 
problem when a = —1, or Oc ^ —1.18 for this parameter set (6 = 1, r = 1, and D = 0.3). The 
solid line follows from our approximate determination of p{x) in the limit of small r. There are no 
adjustable parameters. 
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FIG. 4: Numerically determined bifurcation diagram for Eq. ([T]) with r = 1 (o) defined as the 
point in the (a, b) plane for which a = —1, the exponent of the stationary probability distribution 
function. For reference, we also show the exact bifurcation diagram for the deterministic equation 
^ = (solid line) , and the exact results for the bifurcation threshold of the first moment (x) from 
the linearized equation P] (dashed line). The two points labeled by + lie on the line 6 = 0, and 
are known results for the case of no delay for (from left to right) (x), and for deterministic case 
of ^ = 0. The dotted line is the approximate threshold [Eq. (123]) ]. and the labels (o) are the 
numerically calculated Hopf branch from the probability distribution function of the maximum 
amplitude of the Fourier transform of the trajectories. 
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FIG. 5: Bifurcation threshold Oc from p{x) as a function of noise intensity I? for 6 = 1 and r = 1. 
Time averages used for the determination of p{x) are in t G (300,350), and 1 x 10^ independent 
reahzations have been considered. The hne in the figure is the prediction from our approximate 
determination of the stationary probabihty distribution function. There are no adjustable param- 
eters. 
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FIG. 6: Probability distribution function of the maximum amplitude of the Fourier transform 
Ps{max[X(u;)]) for b = —2, r = 1 and D = 0.3 at the times given, and averaged over 10^ indepen- 
dent realizations. Values of the control parameter shown are: (a) a = —1.1 with a ~ —1.43, (b) 
a = —0.9 with a ~ —0.70, and (c) a = —0.5 with a ~ 0.02, as shown by the solid lines, placed 

above the curves for clarity. The solid lines extend over the range used for estimation of a. The 
distributions in (a) show a clear transient, whereas those in (b) and (c) are stationary. 

19 




FIG. 7: Results of a power law fit to p{x) (o) and to the probability of the maximum amplitude 
of the Fourier transform ps{max[X{u!)]) (o), for small x or small max[X{u!)]. The fitted value 
of the exponent a is shown as a function of the control parameter a. We define the bifurcation 

threshold for the stochastic problem when a = —1, or Oc — —0.999 (o), and Oc — —0.981 (o) for 
this parameter set (6 = —2, r = 1, and D = 0.3). 
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FIG. 8: Ensemble average of x{t — r) given x{t), {xT\x{t)). We have set a = 1,6 = — 0.25,r = 1 
and D = 0.3 The average is computed over 10^ reaUzations for a time interval t G (300,350). The 
dashed line is {xr\x{t)) = (1 + r(a + 6 + D))x{t), the approximate analytic result in the limit of 
small r. 
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